www.gusucode.com > 循环自相关函数工具箱源码程序 > matlab代做 修改 程序循环自相关函数工具箱/cyclostationary_toolbox/frac_delay_cyclic_ac.m
function R=frac_delay_cyclic_ac(x,N,alpha,max_tau) % % FRAC_DELAY_CYCLIC_AC % % calculates the cyclic auto correlation for signal x % at frequency alpha resampling 1/alpha to N samples % max_tau is expressed at original sampling fequency % % R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)} % for k=0 ... N-1 % % USAGE % R=frac_delay_cyclic_ac(x,N,alpha,max_tau) % % calculate cyclic autocorrelation up to max_tau time lags % File: frac_delay_cyclic_ac.m % Last Revised: 5/5/98 % Created: 5/5/98 % Author: Andrew C. McCormick % (C) University of Strathclyde lx=length(x); Rxx=zeros(lx-2*max_tau,2*max_tau+1); for k=-max_tau:max_tau Rxx(:,k+max_tau+1)=[x(max_tau+1:lx-max_tau).*x(k+(max_tau+1:lx-max_tau))]'; end r=zeros(N,2*max_tau+1); w=zeros(N,1); T=1/alpha; for k=1:lx-2*max_tau cycle=floor(k/T); position=(k-cycle*T)/T*N+1; fp=floor(position);cp=ceil(position); if (fp==cp) if (fp==0) w(N)=w(N)+1; r(N,:)=r(N,:)+Rxx(k,:); else w(fp)=w(fp)+1; r(fp,:)=r(fp,:)+Rxx(k,:); end else if (fp==0) w(N)=w(N)+position-fp; r(N,:)=r(N,:)+Rxx(k,:)*(position-fp); else w(fp)=w(fp)+position-fp; r(fp,:)=r(fp,:)+Rxx(k,:)*(position-fp); end if cp>N w(1)=w(1)+cp-position; r(1,:)=r(1,:)+Rxx(k,:)*(cp-position); else w(cp)=w(cp)+cp-position; r(cp,:)=r(cp,:)+Rxx(k,:)*(cp-position); end end end % If insufficient samples to estimate r, fill in blanks with % the neighbouring sample in r. if (sum(w==0)>0) warning('Insufficient samples for desired resolution'); for k=1:N if (w(k)==0) if k>1 r(k,:)=r(k-1,:); end else r(k,:)=r(k,:)/w(k); end end else r=r./(w*ones(1,2*max_tau+1)); end % Take DFT of time varying correlation with appropriate phase change % to compensate for time shift R=zeros(2*max_tau+1,N); for k=-max_tau:max_tau R(k+1+max_tau,:)=exp(-i*pi*((0:N-1)/N)*k).*fft(r(:,k+1+max_tau)'); end